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Effective procedures are presented for the response analysis of 
the SSME turbopumps under transient loading conditions. Of particular 
concern Is the determination of the nonlinear response of the systems 
to rotor Imbalance In presence of bearing clearances. The proposed 
procedures take advantage of the nonlinearities Involved being local- 
ized at only few rotor/housing coupling points. 

The methods Include those based on Integral formulations for the 
Incremental solutions Involving the transition matrices of the rotor 
and housing. Alternatively, a convolutional representation of the 
housing displacements at the coupling points Is proposed which would 
allow performing the transient analysis on a reduced model of the 
housing. The Integral approach Is applied to small d3mamlcal models 
to demonstrate the efficiency of the approach 

For purposes of assessing the numerical Integration results for 
the nonlinear rotor/housing systems, a numerical harmonic balance 
procedure Is developed to enable determining all possible harmonic, 
subharmonic and nonperiodic solutions of the systems. A brief account 
of the Fourier approach Is presented as applied to a two degree of 
freedom rotor-support system. 
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I. 


INTRODUCTION 


The dynamic transient response analysis of the SSME turbopurops is 
essential for further development and prediction of the engine perfor- 
mance under various load levels and required maneuvers* Childs [1] 
conducted numerical analysis for critical speed -transit Ions of the 
HPOTP* A rubbing condition was predicted for that earlier configura- 
tion at the turbine floating -ring seals during shutdown. More recent- 
ly, Childs [2] concluded a series of studies concerning the develop- 
ment of a reliable RPL engine and a description of new problems which 
are being encountered in developing FPL performance* The analysis was 
based on a modal method developed earlier by the same author [3]* 
Childs [2] demonstrated that although a linear transient analysis 
remains an efficient procedure for general characterization of the 
turbopump’s rotor -dynami cs , nonlinear analysis is essential* An 
important case considered by Childs is that of the effect of the 
radial clearances provided at the outer races of the bearings* The 
results of the nonlinear analysis showed a significant reduction in 
the subsynchronous rotor motion* More significantly, the bearing 
clearances can drop the peak -vibration running speed into the 
operating range where the synchronous whirling loads might pose a 
serious threat to bearing life* 

For large complex rotor systems, such as the SSME turbopumps, 
various modeling and analysis techniques vary in their ability to 
accurately describe the systems’ behavior* This ability mainly 



depends on the configuration of the systems analyzed and on the for- 
cing conditions. In rotating components, this also Involves whether 
the rotating speed is constant or varying with time. A hybrid repre- 
sentation by various types of coordinates and formulation for the 
various components of the systems may prove valuable. 

Different procedures have been utilized by analysts to determine 
the transient response of large order rotor systems. The procedures 
can be recognized as falling under one of two basic approaches. Those 
using physical or modal coordinates of the complete system and those 
using the coordinates of the individual components of the system. The 
methods also differ in the numerical integration methods selected for 
the analysis. 

Rouch and Kao [4] employed Guyan (static) reduction method to 
arrive at a reduced size model in terms of the remaining physical 
coordinates. Accuracy of the results could be expected to be accept- 
able since the rotor is basically a train of mass -stiffness subsys- 
tems. Nordmann [5] attempted to minimize the Inaccuracy of static 
condensation by applying the static reduction technique to an arbi- 
trarily substructured rotor system and then assemble the reduced sub- 
structures to form a reduced system. The procedures is very laborious 
and no guarantees of accuracy are apparent. 

Childs [6] utilized free -interface modes of the various system 
components to represent the assembled SSME turbopumps. The method, 
using fourth order Runge-Kutta integration, does not provide for 
accommodating accurate modal representation of the large housing model 
while maintaining a small size for the model. Nelson e^ al^. [7] on 


the other hand used fixed-interface complex component modes to assem- 
ble a reduced size model. For systems with large number of coupling 
points among the components, the approach suffers the problem of 
introducing higher frequencies resulting from excessive number of 
constraints Imposed at the coupling (or boundary) points. In a tran- 
sient analysis, this will necessarily result In much smaller time 
Increments and consequently, will lead to excessive computational time 
and larger round-off errors. 

For nonlinear large rotor systems, only a few analysts have pre- 
sented techniques for the general transient analysis of such systems. 
Adams [8] used a normal mode representation for the rotor In terms of 
its undamped, free s)TDmetrlc modes and treated gyroscopic and non- 
linear terms as pseudo -external loads. The method presented by Childs 
[6] makes use of a similar procedure to couple the rotor to Its flex- 
ible housing. Nelson ei^ al . [7] developed a general computer code for 
the transient analysis of large rotor systems. The user may utilize 
time-step Integration In the constrained-rotor (fixed-interface) modal 
space. Again, all connection points. Including those at the non- 
linearltles must be constrained, leading to the same shortcomings 
described previously. 

None of the above studies have adequately addressed the problems 
of attempting to use reduced size, accurate models and the associated 
efficiency of computation. A judicious selection of the system con- 
figurations and the numerical methods Is essential for achieving the 
required efficiency while maintaining acceptable accuracy. 



In this study, integration methods based on the transition 
matrices [9] of the separate rotor and housing are used to efficiently 
determine the nonlinear transient coupled system response under imbal- 
ance forces and in presence of bearing clearances* Pragenau [10] 
utilized the transition matrix for integration, stating that it offers 
the simplicity of the Euler method without requiring small steps. 
Pragenau maintained that for constant subsystems, the stability and 
accuracy of the method are acquired through the closed form solution 
of the transition matrix. 

Alternatively, a discretized Duhamel [9] (convolution) integral 
method can be used to an advantage to represent the response at the 
housing coupling points to the rotor. Kubomura [11] used a convolu- 
tion based method to achieve dynamic condensation of a substructure to 
its coupling points to other structures. Convolution was also used in 
[12,13] to reduce system coordinates to that of the nonlinearities. 

Along with the response analysis in presence of bearing clear- 
ances under Imbalance forces, a numerical harmonic balance method is 
developed toward verifying possible steady state synchronous and sub- 
synchronous rotor response. This is essential to ensure that no 
possible potentially damaging solution is missed due to unfortunate 
selections of particular initial conditions. The harmonic method 
locates all possible periodic solutions. The method is briefly 
presented here from references [14] and [15] as applied to a modified 


Jeffcott model 



II. THE TURBOPUMP MODEL 


Housing: 

The modal equations of the housing In the X-Z, Y-Z planes In 

terms of a truncated set of Its modal coordinates {q„}, normalized 

H 

with respect to the mass matrix, can be written as [16,17] 

{q^} + ['2 Cj, {qjj} = (D 

where ^ and A are the modal damping coefficients and natural fre- 
H n 

quencles, respectively, is the vector of coupling forces to the 
rotor, including the balance piston axial force. The axial force Is a 
function of the axial displacement and velocity as well as the spin- 
ning speed. The coupled physical displacements in the X and Y direc- 
tions at the coupling points to the rotor are given by 

{qg} (2) 

where [Ag(j] is the normalized modal vectors, with respect to the 
mass matrix, associated with the coupling points. 

Rotor: 

As In the case of housing, the symmetric -rot or equations of 
motion may be written In terms of a truncated set of modal coordinates 
{qg} of the free-free nonspinning, undamped rotor as 
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• where and are nondiagonal matrices which are speed, ([i, depen- 
dent, and Is the Imbalance forces which Is, In general, functions 

• • • 

of (|; and (|;. The vector represents the Interaction forces with the 

• housing such as the bearing, seal, side loads (which are functions of 

, etc# the physical displacements at these coupling points are 
given by 


K AfT 

{w } = {„^} = [\c 


.?T- ^‘’r^ 

^RC 


(4) 


Interaction Forces 


The coupling forces between the housing and rotor can be 
expressed as (see Fig. 1) [17] 

- {Fjj} = [C] {Wjj} + [K] {Wjj} - [C] {Wg} - [K] {W^} 


and 


{Fr} = {Fj^} 


(5) 


( 6 ) 


where [C] and [K] Involve direct and cross coupled stiffness and damp- 
ing forces as well as spinning velocity dependent coefficients. The 
coefficients of the bearing forces allow for presence of clearances 
(deadbands) . 







8 


III. TRANSITION MATRIX FORMULATIONS 

The rotor's modal equations of motion are written In the first 
order form 


<®R> - >“r) <“r> * ‘"r* 


(7) 


where 


<”r> 


<,:> > <“r> = -“r, 


( 8 ) 


and 




(9) 


Similar first order equations can be written for the housing. 

The solution of equations (7) with constant rotating speed and 
the analogous ones for the housing In terms of the associated transi- 
tion matrices of the rotor and housing takes the form [9] 


(U(t)} = {U(0)} + q/ e^“^^‘^~”'^{F(T)} dt (10) 


This representation can be cast In discretized form as 




( 11 ) 


where {Uj} = (U(tj)} 


and T = tj+j - tj 


( 12 ) 



The force vector In eq. (11) can be either treated as (1) con- 
stant within a small time Increment T or (il) a linear function of 
time. Clearly, the linear representation would result In more accuracy 
for a given Increment T. An assumption of a step load (constant load 
within an Increment) allows equation (11) to be written in the simple 
form 


(Ui+i) = [$(T)] {Uj} + ([f(T)J - ['!.]) [a] ^ (Fi) 

where t4^(T)] = , {Fj} = {F(ti)} 


^[a]T 

e 


CO 


1 

n=0 


n! 


[a]" 


(13) 

(14) 

(15) 


so that ($(T) -I) [a]“^ + I 7 ^;^, [al"1 (16) 

n=l 

Both expressions in (15) and (16) converge very rapidly for small 
Increment, T, and when used as above (constant [a]), need only be 
calculated once for the entire time response history of the coupled 
rotor/housing system. 

A more efficient algorithm can be constructed by representing the 
coupling forces as linear functions of time within each increment T, 
or 

{F(t)} = {Fj} + — Y~^ (<fi+l} - {Fi>) (17) 


Using this representation, and after some manipulation, the solution 
(11) can be shown to take the form [10] 
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{Ui+i} = [$]{Ui> + [[#] - ['Ijl [a]"^ ({Fi> + [a] 


1 1 + 

^ r„ •> . r i“l 1 + 1 1 •> 

^ 7 


- [a]'^ ({Fi+i} - {Fj}) 


(18) 


An incremental solution using equation (18) for the rotor and a simi- 
lar matrix equation for the housing, along with equations (5) and (6) 
in first order form, can be employed to construct the response time 
history of the turbopump considered. The particular recurrence proce- 
dure used will depend on the type of transient response sought as well 
as on the accuracy versus computational time tradeoffs. Some discus- 
sion of the computational procedure is presented in a later section. 

Use of The Duhamel Integral 

An efficient incremental representation of the housing displace- 
ments can be achieved using the convolution Integral. As with the 
transition matrix formulation, the coupling forces are treated as 
external loads on the housing and may be assumed linear in time for 
every Increment T. The response at any time t is given for zero 
initial conditions, using equation (1), as 

t -E wCt-x) 

{Wjj}= / ['-i e sin w^Ct-x) J [A^^] {Fy( t) M t (19) 

o (i 


or 


(Wjj(t)} = 


w ^ ^ (t-x) 

, 1 , 

j=l o “5 

sin 0 ^ (t~x) {F( x) } d X 


( 20 ) 



Where (A Is a reduced column matrix of elements of the jth 

eigenvector associated with the coupling points of the housing and 
rotor, is the damped natural frequency of the nth mode and "tr” 
stands for a transpose. 

To enhance the accuracy of the modal representation of the hous-- 
Ing, the deleted higher housing modes can be approximately accounted 
for through the residual flexibility corresponding to these modes, or 
at any time t, 

<"h' - ’ ‘«H> '■'h* ■ »» 

where [G^l Is the residual flexibility matrix corresponding to the 
coupling points on the housing's nth mode of the uncoupled housing. 

For the purpose of demonstrating the method In simpler terms, 
consider the generalized housing coordinate of the undamped nth mode 


«!n(t) 


1 




t 

/ sin Wj,(t-T) Pn(''^) di 
o 


( 22 ) 


where Pn('^) Is the unknown modal coupling force. If the force Is 
assumed linear In time within each Increment, the generalized dis- 
placement In equation (22) due to a coupling force applied between 
tj and tj+j is given by 

! IPn(ti) + -=r- (Pn(ti+i)-Pn(ti)J}sln a)„(t-T)dT (23) 
“ ti 

After some manipulation, the Integration in equation (23) can be 
written in closed form In terms of the unknowns P^'s as 
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qn(t) = {cos ojnt (Pn(ti + l) cos ^ntl+l ~ PnCtj) cos 


Wr 


V'iti) - V'l' 


03 , 


n 


(sin - sin Wjjti)) 


+ sin Ujjt (Pn(ti+i) sin - PnCtj) sin 


P (t ) - P (t ) 

+ (cos Wjjti+i ~ cos ^ntl))} 


(24) 


or 


qjj(t) = — ^ cos tOjjt + sin Wj^t) 




(25) 


The total response at time t due to the contribution of the coupling 
forces from zero to time t (= N.T.) Is 


qn(t) = 


1 


Wr- 


[cos Wjj NT I 
1=0 


.( 1 ) 


N . 

sin coj, NT I 


(26) 


1=0 


The generalized velocity can be expressed as 


qn(t) = ^ [cos WnN.T I sin N.T I (27) 


1=0 


1=0 


The generalized acceleration can be obtained from 

qn(t) = -^n qn(t) + Pn(t) (28) 

The generalized displacements and velocities of the rotor In Incremen- 
tal form are given previously In the form of equation (18). Equations 


(18), (27), (5) and (6) can be readily employed to calculate the 

response time history for the coupled rotor/housing system* 


Computational Procedures and Results 


Various alternative procedures have been explored for the 
efficient implementation of the integral formulations based on the 
transition matrices and the convolution integral. Although studies 
are continuing to fully make use of the computational advantages 
offered by these formulations, the following appears to be the most 
attractive to date. 

For both the transition and convolution methods, a linearized 
representation of the coupling forces (including those at the bearing 
deadband) are utilized. Equations of the form of (18), (26) and (27) 
are used along with the coupling forces relations given by (5) and 
(6). If the rotor system is linear , the forces are expressed in terms 
of the displacements at the coupling points. The system is then 
represented by a simultaneous system of equations involving the hous- 
ing and rotor coupling coordinates so that at each increment of time 



> = (f(W^ 



(29) 


where f stands for a function. 


Another alternative which proves to be very effective specially 
In presence of the bearing deadband nonlinearity is by again using a 
lineared force variation within each time Increment, but including an 
Iterative procedure at each time step. The iteration work as follows. 
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where F stands for coupling forces: Predict the response at the 

coupling points at tj 4 j from those at tj 


“ii.l • ‘h 


Ri+1 






Hi 

R 


Imbalance forces) 


(30) 


Then predict the forces at tj+j from 


F' 

Hl+l 

Hl+l 



(W* 

Hl+l 


(W 


H 


1+1 



) 


W' , imbalance forces) 
Hi+l 


(31) 


Correct the displacements at tj+i as 


W 

Hl+l 





w 

Ri+1 


R 




F’ 

Hi+i 


) 


F' , Imbalance forces) 

Hl+l 


(32) 


Re-iterate using equations (31) and (32). If the maximum differ- 
ence between the displacements in two consecutive iterations Is 
smaller than a specified tolerance, the Iteration Is stopped and the 
final values of displacements at t^^j taken as 


W s U 

Hi+i %+i 

K “ Wp 

Rl+l Hj+i 


(33) 


The transition matrix method was applied to a system consisting 
of damped two subsystems of three and two degrees of freedom Interact- 
ing through a clearance* The results are presented In the Appendix, 
in which comparisons to the results obtained using an iterative Runge- 
Kutta method are made* The results show the transition matrix method 



to be more efficient than that of the numerical integration procedure 
for a given time increment* However, more work is needed to devise an 
alternative iterative technique which would expand the range of 
convergence of the method* 

For a subsystem to which the form of the convolution integral 
(19) is applicable, a solution procedure similar to that involving 
equations (30) to (33) can be applied using equation (26)* 

Alternatively, a one step procedure can be utilized if the coup- 
ling forces are assumed to be constant, rather than varying linearly, 
within each time step* In that case, the expression for qn(^) is 
obtained by setting 

The efficiency of this formulation can be greatly enhanced by util- 
izing a Taylor series expansion for the displacements at time 
terms of the displacements and their derivatives at t^» or 

"tl*! ■ “t, * “ti <^1-1 - 

2 “cj "l+l ■ 

O DO 

where and are given by appropriate expressions of the forms 

given by equations (27) and (28). The displacements (and velocities) 
as In equation (34) are used to calculate the coupling force In each 
time step. 

An approach combining the Duhamel representation for the housing 
of a rotor system and a transition matrix algorithm for the rotor 



could be utilized to an advantage in determining the transient 
response of the system. The approach is yet to be applied to rotor- 
housing systems. 


IV. NUMERICAL HARMONIC BALANCE METHOD 


A method is developed [15] In order to check all possible nonlin- 
ear responses of simple rotor systems under periodic Imbalance loads 
In presence of bearing clearances, rubbing, seal forces, side forces 
and others. The method requires the spinning speed of the rotor to be 
constant so that the model is represented by a nonhomogeneous system 
of nonlinear equations of constant coefficients. The developed method 
is a modification of that due to Yamauchl [18]. 

The early work in nonlinear rotordynamics by Yamamoto [20] intro- 
duces a nonlinearity to the Jeffcott equation by Including the effect 
of bearing clearances. The Van der Pol self -sustained vibration 
theory was used for the analysis under the assumption of small clear- 
ances. Due to the symmetry of the gaps, only the harmonic response 
was considered. Recently, Childs [21] presented an explanation for 
the subharmonic response of rotors in presence of bearing clearances 
and a side load, using perturbation techniques under the assumption of 
small nonlinearity. The potential destructive vibration of the LOX 
pump in presence of rub in the SSME (Space Shuttle Main Engine) has 
been considered by Beatty [22]. 

Glease and Bukley [23] used a direct numerical integration method 
to analyze the rotor response, including seal and nonlinear bearing 
forces. However, the results only show the transient response before 
a steady -state is reached. Since the damping is usually small in 
rotor systems, a steady-state response is often established after a 
relatively long period of integration time. For the steady-state 
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response, Yamauchi [18] developed a numerical harmonic balance method 
using the FFT (Fast Fourier Transformation) algorithm for nonlinear 
multiple degrees of freedom rotor systems, in which transfer matrix 
formulations were used to describe the system. Salto [19] calculated 
the nonlinear unbalance response of horizontal Jeffcott rotors. Both 
studies of Yamauchi and Salto were concerned only with the harmonic 
response. No analysis for the Important case of subharmonic rotor 

response was attempted. 

The procedure developed in [15], and presented here, is based on 
a modified numerical harmonic balance method using discrete Fourier 
transformation and its inverse. This transformation, rather than an 
FFT procedure, was utilized in order to reduce computational time and 
errors. This is achieved by calculating the complex exponential 

values at the beginning of computation and then store them in active 
memory for subsequent calculations. Subharmonic and superharmonic 
responses are also accounted for by the method. An account of the 
method and some results obtained for the response of the modified 
Jeffcott rotor system selected for the analysis is outlined in what 
follows m 

The Modified Jeffcott Model 

The rotor model is depicted in Fig. 2. In the figure, ^ is the 
shaft speed which is assumed constant, and ^ is the angle of the disk 
rotation with respect to the inertial coordinate system y - z and is 
referred to as the whirl angle. The displacement of the rotor center 
from the origin of the inertial references frame is denoted by r. The 
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eccentricity, e , denotes the displacement of the mass center from the 
geometric center of the shaft. 

The forces considered in the formulation of the equation of 
motion Include bearing forces, seal forces, an imbalance force and a 
side force. The bearing forces can be assumed as plecewlse-linear and 
occur only when the displacement of the shaft is greater than the 
clearance (or deadband) between bearing outer race and support. Fric- 
tion due to rubbing between shaft and bearing support during contacts 
is also considered. The bearing forces, F^, can be expressed as 
follows. 

Kb(r-6uj-) + pK^jUx X (r-6ur) ; jr| 2 

Fb - (35) 

0 ; r < 6 

-► -► 

in which and u^ are unit vectors in the r and x direction, respec- 
tively, Kb is the bearing stiffness, 6 is the deadband and p is the 
coefficient of friction between shaft and bearing. The "X" stands for 
a vector cross product. The assumed form representing the seal forces 
is given by 

Where Kg is the seal direct stiffness, Cg is the seal damping coeffi- 
cient, and Qg is the seal cross-coupling stiffness. The side force is 
assumed to be due to gravity. The gravity loading provides a side 
force In the negative direction. 

The force equilibrium equations for the rotor can be written, 
with respect to the y - z inertia reference coordinates, as follows: 
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My 


Cgy 


Key 


+ QgZ + Kb (y - y 


/ y + 


:) 

r 


- nKb (z - z 


y 2 2 

r y + z 


) I = Mew2 coswt - Mg 


(37-a) 


• • • r 

Mz + CgZ + KgZ - Qgy + K^, (z - z — ) 

O O 


^2^2 
/ y + Z 


+ ^Kb (y - y 


y 2 ^ 2 

/ y +2 


) = Mew^slncot, 


(37-b) 


• in which M Is the mass of the disk, g Is the acceleration of gravity, 

and I I denotes that the nonlinear bearing forces occur only when 

/y2 + z2 > ; otherwise, they are zero* In the above equations, the 

® rubbing between stator and blade tip Instead of bearing rubbing (which 

Is usually small) can be calculated by changing Kb to that of the 
stator stiffness. 

® The equations of motions, (37-a) and (37-b), can be put In 

nondlmentlonallzed form using the deadband, 6, as the reference 
displacement and the natural frequency of the associated linear prob- 

• lems (6 = 0) as the reference frequency. This natural frequency Is 
not that of the nonlinear system, although due to small clearance the 
periods of maintained contact during rotation may lead to a frequency 

• close to that of the associated linear frequency. Let 


Kq = Kg + Kb 


2 

0) 

n 


M 


Q 

“n 






( 38 ) 
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Y 


= Z 
6 



R = 






d) = 


_s_ 

0)^6 

n 


v0 = U)t 

N 

where 0 is a normalized time, and v is the subharmonic order. Using 
the parameters in equation (38) with equations (37*a) and (37-b) leads 
to 


r* 



Y’ 


+ 


CtV V 9 V 

— -Y + y —Z + G(0) - HF(0) = ev^cos(v0) + 4^ 


(39-a) 


- 2 2 

Z" + + -^Z - Y -^Y + F(0) + |iG(0) 


eV^ sln(v0) 


(39-b) 


where 


G(6) 


«v2 

P^Y 


(1 



(39-c) 


v2 , 1 , 

F(0) - P^Z|(1 --|)| (39-d) 

in which a prime denotes differentiation with respect to the normal- 
Ized time 6 and the nonlinear restoring forces are established when- 
ever R is greater than unity. 

Method of Analysis 

The steady periodic solution of equations (39), including any 
subharmonic, superharmonic or harmonic vibration can be written as 


Fourier series, or 
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# 




N 


Y(0) = ayO + 2 1 

n=l 

(ayn cos 

n0 

■“ n0 ) 

(40-a) 

N 

z(e) =820+2 1 

n“l 

(a 2 n cos 

ne 

- b2nSln ne) 

(40-b) 


where N is the number of harmonics to be taken into account in the 
final solution under the assumption of small frequency bandwidth. The 
multiple of two and the minus sign associated with the bn are adopted 
to facilitate accommodating complex Fourier coefficients when discrete 
time data is used. In the same say, the time series representation of 
the nonlinear bearing forces, occurring only when j R | > 1, can be 
expressed as Fourier series of the form: 

N 

G(0) = CyO +2 J (cyn cos n6 - dynsln n0) (41-b) 

n=l 

N 

F(0) = C 20 ■•■2 'I (C 2 n cos n0 - d 2 nSln n0) (41-b) 

n®l 

Since these nonlinear bearing forces are a consequence of the exist- 
ence of a deadband in the system, the Fourier coefficients of the non- 
linear restoring forces are functions of the Fourier coefficients of 
the steady periodic solution sought. Substituting equations (40) and 
(41) into equations (39-a) and (39-b) and applying a harmonic balance 
procedure yields 4N + 2 nonlinear simultaneous equations involving 8N 
+ 4 unknowns. The harmonic balance procedure gives for the constant 


term: 



24 


e 


2 2 2 
gy . j. V .V 

~Y- ayO + CyO + Y ~2 3zo " “ ^P~T 

C Q Q 


(42-a) 


2 2 

2 ®zO ^zO ~ Y 2 ^yo ^ ^^^y0 * ® 
C Q 


(42-b) 


the cosine terms: 


2 2Cv 2 

-n ayn - n byn + o~ ayn + Cyn + ^zn ~ * O.Smev (42-c) 

Q Q 


2 2^ 
1 a,„ - n 


®zn “ ^ ^zn ■*■ tt~2 "*■ ^zn 

Q 


C 


2 ®yn '*’ M'^yn ~ ® (42~d) 


and the sine terms: 


2 

” byn “ n q ayn “ byn - dy„ - bzn + |idzn 


V 


= 0 


(42-e) 


2, 2Cv 

n bzn - ^ Q ^zn 


®~2 “ ^zn + r ~2 byn ~ ^^dy„ = O.Smev (42-f) 


in which 


ID - 0 ; n ^ V 

m = 1 ; n *= V 


Another set of 4N + 2 equations can be found from the relation- 
ships between the Fourier coefficients of the nonlinear bearing forces 
and those of the steady-estate periodic solution. These relationships 
can be found by determining Y(0) and Z(9) from the Fourier coeffi- 
cients of the steady-state periodic solution using the discrete 
inverse Fourier transformation, 
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N-1 

dp “ Resl( J 
k=0 


l(2nkr/N)' 


where Dj^ stands for the Fourier coefficient, and dp, for discrete dis- 
placement. Using equations (35), discrete time series of the nonlin- 
ear bearing force can be found. The inverse discrete Fourier trans- 
formation of the time series solution for the nonlinear bearing forces 
will yield the Fourier coefficients of the nonlinear bearing force as 


Bk - (V 

r*0 


The resulting nonlinear simultaneous algebraic equations can be 
handled using a Newton -Raphson method. The Newton-Raphson method uses 
an incremental procedure in determining the values for the next itera- 
tion as follows. Let 


S « S + As, 


where S represents the Fourier coefficients of the steady-state solu- 
tion, and the superscript ”0’* denotes current state while As stands 
for the increments of the coefficients during one step of iteration. 
For example, 

^yn “ ^yn 

Similarly, the Fourier coefficients of the nonlinear restoring 
force can also be expressed as follows: 


B + B + Ab, 
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For example, 


^yn ‘^yn ^*^yn 


In which the Increment AB can be expressed using the total derivatives 
with respect to all the Fourier coefficients of the steady-state solu- 
tion, S®, so that 


For example, 


- I ^ iSi. 

n=l 


N 5c 6c 


(A7) 


N 6c 


Acyn = I (-^ Aayn + -j^Aa ) + j (-^ Ah 

n=0 ° yn °®zn zn dby^ 


6c 


zn 


where the partial derivatives are to be calculated at the current 
state value. 

In order to account for the large nonlinearity in the system, the 

increments, ASn, must be chosen small as compared to whenever 

numerical differentiation is performed. The numerical differentiation 

is performed using forward differentiation. For example, to obtain 
6c 


the derivative 


JUL 


“yn 


, a new Cyn is calculated at ayn = ayn + Aayn from 


»yn 


*yn 


the nonlinear relations of equations (39-c) and (39-d). The numerical 
differentiation is then given by 

0 


6c 

5a 


c - c 
yn yn 


‘yn 


Aa 


yn 
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The numerical differentiation procedure can be done using forward 
differentiation as follows: 


1 . 

2 . 

3. 

4. 

5. 

6 . 


3q Is changed by adding the preset small Aao» other 

Fourier coefficients, a^ and bn» are not changed. 

The newly assumed discrete displacements, Y(G) and Z(0), are 
generated using a discrete Inverse Fourier transform. 
Discrete nonlinear forces, G(6) and F(0), are calculated 
using equations (39). 

The Fourier coefficients of the discrete nonlinear force, 
co» Cn» 3nd dn> are calculated using Discrete Fourier Trans- 
formation. 

The differential values of -r — , - — -, and are calculated 

Sao dao dao 

The calculations In step 1 to 5 are repeated for each a± and 


bf. 


Substituting equations (45) and (46) Into equations (42), the 
following Incremental form will result. 


2 2 

o ^ayo + ^CyO y o ^®zo ~ 
Q Q 


2 2 - 
V ov 0 


0 


t 2 2 *^y® ~ ^ 2 ^ ^^^zo 


Q Q 


(48-a) 


2 2 
(XV V 

— J AazO + ACgO - Aayo + pACyQ 

Q 


7 7 

av 0 0 ^ V 

— 5 azo - CzO + Y “o 


»yo 


0 

^CyO, 


(48-b) 
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for the cosine terms; 


2 2Cv 

n Aayn “ n Q ^yn ■*■ ^®zn “ ^^^zn 


0.5 mev^ + n^Syn + " byn - 3 yn - Cyn - ajj, + jJcJn 


2 

”^^zn " n “o“ ^zn + ^zn ~ ^vn + MAc 


V V 

-Q- “‘'zn ^zn + Aczn ~ Y-J “«yn 

£c Q 


yn 


2 0 

= n azn + n 


2 Cv ,0 

IT ^zn 


.v2 0 




*zn 


^zn ■*■ 


0 

«yn 


Y-^ avn ~ Me 


yn 


and for the sine terms: 


^^^^yn “ a Q ^ayn “ ® 2 ^^yn ~ ^<^yn ~ Y— j ^^zn ■•■ MAd 


Q" 


zn 


2. 0 ^ 2Cv 0 . . 0 jO ,0 ,0 

- n byn + n -jj- ayn + a-^ byn + dy^ + y— ~ Md^j, 


2 a 2 Cv , 

n Ab^n - n Aazn " ^l>zn ~ ^<^zn ■*■ Y-^ Aby„ - JxAdyn 


0.5 mEV - 


2,0 ^ 2Cv 0 ^ ^0 ,0 

n t>2n + n ^zn + “-, bzn + dzn 


Y^ byn + Mdyn 


In which 


m “ 0 
m = 1 


n ^ V 
n = V 


Equations ( 48 -a) to ( 48 -f) can be put In the matrix form 


( 48 -c) 


( 48 -d) 


( 48 -e) 


( 48 -f) 


[K] {Ax> = (w) 


( 49 ) 



Where [K] corresponds to the Jacobian matrix whose elements are evalu- 
ated at every step. 

This Iteration Is continued until all the components of the 
correction vector (w) become sufficiently close to zero. The calcula- 
tion procedure presented here Is basically an Iterative procedure. By 
setting the Initial values to zero, a linear solution Is calculated In 
the first Iteration. The linear solution Is modified as the Iteration 
progresses. For the next frequency ratio, the final solution of the 
previous step frequency ratio Is used to facilitate the convergence of 
the Iteration. 

Numerical Results and Discussion 

A comparison between results using direct numerical Integration 
and those of the FFT method Is shown In Fig. 3. The numerical Inte- 
gration Is carried out using the central difference method with the 
nondlmenslonal equations (39-a) and (39-b). To find the steady-state 
solution within a reasonable duration, the damping factor, C, Is 
chosen to be relatively large. For the application of the FFT method 
In a numerical procedure, four harmonic coefficients are taken Into 
account for each of the displacements, Y and Z. Also, sixteen dis- 
crete data points are chosen In one period to avoid the aliasing 
effects. The mean displacements, Y, Z, and R, are defined as half of 
the sum of the maximum displacement and the minimum displacement In a 
steady-state response. As shown In Fig. 3, the frequency response 
results using the FFT method demonstrate good agreement with those of 
the direct numerical Integration method. 







In Fig. 4, the nonlinear response of the system is shown in pres- 
ence of the side force due to gravity. To ensure the occurrence of 
subharmonic response, the damping factor, C, was set to zero* As 
depicted in Fig. 4, a second-order superharmonic resonance occurs near 
the frequency ratio of 0*5. A second-order subharmonic resonance 
occurs below the frequency ratio of two, while a third -subharmonic 
resonance appears below the frequency ratio of three. These trends 
have already been demonstrated experimentally by Bently [26] and Erich 
[27]. The nonlinear steady-state analysis also reveals the co-exist- 
ence of harmonic response with the subharmonic response. The occur- 
rence of a specific type of response depends upon the particular 
initial conditions (due to proximity to a corresponding domain of 
attraction) . 

A comparison between the nonlinear response and harmonic motions 
in the y-z plane can be made by examining the results presented in 
Fig. 5, 6, and 7. The shaft center traces noncrossing paths in a 
harmonic response case, or otherwise the paths will be of more compli- 
cated shapes and larger radii (a less desirable behavior in rotating 
machinery). The existence of second -order subharmonic resonances is 
examined for various values of the nondlmenslonallzed side force 
factor, 0 * in Fig. 8. The figure shows that subharmonic vlbra- 

tion would not exist for either zero or relatively large values of •t’. 
This Is because only synmetrlc motion Is maintained for these particu- 
lar values of <t>. Fig. 8 also shows that a smaller value results In 
a subharmonic response within a broader frequency range. The choice 
of a smaller clearance or a softer bearing stiffness can, therefore. 
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e = 1.5, <s> = 1.5, 
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a = 0.01, P = 1., Y = 0, C = 0, 

P = 0. 


33 




: Second-order superharmonics 






reduce the possibility of a damaging resonance. The effect of eccen- 
tricity can be deduced from the examination of Fig. 9. A larger 
eccentricity Influences the horizontal motion (z-dlrectlon) more than 
it does the vertical motion. 

Other results concerning the effect of rubbing due to friction, 
damping, and cross -coupling stiffness can be found in reference [15]. 


frequency ratio, ft 


Figure 9. Effect of eccentricity on nonlinear rotor dynamic 
response; o = 0.01, P = 1., Y = 0, C = 0.1, 

<t> = 1.5, 11=0 


V 


USE OF COMPONENT MODE METHODS 


A modified fixed interface component mode procedure of Glasgow 
and Nelson [24] can be applied to certain nonlinear rotor systems with 
flexible housing. The nonlinearity is assumed to arise from the exis- 
tence of clearances at the bearings. The coupling between a rotor and 
its housing occurs through bearing, seal, fluid and other interaction 
forces. 

The modified analysis procedure can be carried out as follows . 
The rotor and housing are coupled at their points of interaction 
except at the bearings with deadband clearances. The eigen-parameters 
of the coupled system are obtained with the rotor and housing fixed at 
the locations of the deadbands. The coupled system is then repre- 
sented by a truncated set of these modes plus a static constrained 
mode [24] corresponding to the degrees of freedom at the location of 
the bearing clearances. 

The fixed -Interface component method was applied to two simple 
multi -degree of freedom subsystems and the results showed it to be 
more accurate than the original Nelson ^s method for the same total 
combined number of dynamic and static modes of the system. This is so 
since the method as modified above allows including more numbers of 
dynamic coordinates of the components in the analysis. The modified 
method becomes more efficient than the original method in cases where 
interaction forces occur at relatively large numbers of locations 
between the rotor and housing. The modified approach would allow 
first coupling the housing and rotor using more numbers of modes than 



currently being used, then reduce the number of modes of the resulting 
coupled system. 
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VI. 

1 . 

2 . 


3. 


4. 


DISCUSSION AND RECOMMENDATIONS 

Hybrid component representation and numerical incremental proce- 
dures for the transient response analysis of complex rotor 
systems can lead to more efficient methods. 

The explicit integration methods based on the transition matrices 
and convolution can be very effective in determining the transi- 
ent response of large flexible rotor/housing systems such as the 
SSME turbopumps. The methods are particularly efficient in cases 
concerning constant spinning rotor speeds and in presence of 
bearing deadband clearances. More work is, however, still needed 
to exploit and further develop the methods to their fullest 
potential. 

A modified fixed-interference component mode method could be used 
to construct a reduced size rotor/housing system which is more 
accurate than that of the original method of reference [24]. The 
modification concerns the use of smaller number of connection 
points as the fixed Interfaces of the system. Similarly, the 
hybrid coupling method of McNeal [25] could be extended for 
application to rotor systems. 

A numerical harmonic balance method using discrete Fourier Trans- 
formation is developed and applied to a modified Jeffcott model 
including bearing clearances, seal cross coupling forces, a side 
force and friction due to rubbing. The method can be used to 
determine all possible steady state solutions for the rotor. The 
method can be extended to larger rotor systems, taking advantage 



of the nonlinearities Involved being localized. Application of 
the method will ensure that no potentially damaging periodic non- 
linear response of a given rotor will be missed by solely depend- 
ing on numerical integration methods. Arbitrarily selected 
initial conditions may not necessarily lead to a possible peri- 
odic solution using Integration techniques. The determination of 
the domains of attraction with multiple solutions and the exten- 
sion of the algorithm to multi -degree of freedom rotor systems 
are currently under consideration. 
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APPENDIX 


The Transient Response of a Test Model 

Two single mass spring subsystems Interacting through a gap are 
utilized to form the test model (see Fig. A.l). Modal representation 
of the subsystems Is utilized. The transient analysis Is carried out 
using both the Integral transition matrix method and an Iterative 
Runge-Kutta-Verner Integration scheme. A comparison of the computa- 
tion time of both methods Is shown In Table A.l. The comparison shows 
the transition matrix procedure to be more efficient. 
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TABLE A.l. Comparison of the Solution Time Calculated by Direct 
• Integration and the Transition Matrix Methods 


I. USING ALL THE SUBSYSTEMS MODES 


Time Step 
Size 

Numbe r 
of Steps 

CPU Computation Time 

Runge -Kutta-Verner 

Transition Matrix 

1 X 10"^ 

45000 

15 min. 26.7 sec. 

4 min. 20.5 sec. 

3 X 10"^ 

15000 

10 15 

wamam 

1 X 10"'^ 

4500 

8 46.3 

mmem 

2 X 10‘3 

2250 

8 27 

diverge 

3 X 10-3 : 

1500 

9 1.2 

diverge 


II. USING ONE MODE OF SUBSYSTEM B AND TWO MODES OF A 

1 X 10"^ 

45000 

8 min. 59.8 sec. 

5 min. 22.4 sec* 

3 X 10-^ 

15000 

3 16.6 


1 X 10“^ 

4500 

3 2.8 

1 43.6 

2 X 10~3 

2250 

3 12.3 

diverge 
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Figure A.l. The test model. 








